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Abstract. We investigate transport through molecular wires whose energy levels are affected 
by environmental fluctuations. We assume that the relevant fluctuations are so slow that they, 
within a tight-binding description, can be described by disordered, Gaussian distributed onsite 
energies. For long wires, we find that the corresponding current distribution can be rather broad 
even for a small energy variance. Moreover, we analyse with a Floquet master equation the 
interplay of laser excitations and static disorder. Then the disorder leads to spatial asymmetries 
such that the laser diving can induce a ratchet current. 
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1. Introduction 

Chemical adsorption of sulfur atoms on gold surface allows a stable bond between gold 
tips and thiol groups of molecules. This has been exploited for measuring the conductance 
and the current- voltage characteristics of gold-molecule-gold junctions [1-6]. Repeated 
measurements even with the same sample, however, reveal small but noticeable differences 
which possibly stem from environmental fluctuations that impact upon the effective molecule 
parameters. Moreover, the particular form of the gold tip can have a significant influence on 
the transport properties [7]. 

A present line of experimental research is the measurement of molecular conductance 
when the electrons are excited by electromagnetic waves. There one expects various 
phenomena ranging from photo-assisted transport [8, 9] to ratchet or non-adiabatic pump 
effects, i.e. the induction of dc currents by ac fields even in the absence of any voltage 
bias [10, 11]. Moreover, it has been predicted that properly taylored laser pulses can give 
rise to short current pulses [12-15]. Since a dc current flows into one particular direction, a 
ratchet effect can occur only in "sufficiently asymmetric" systems [9]. In that respect, a static 
disorder is sufficient to break the reflection symmetry of an individual realization and, thus, 
may support a ratchet effect. 

The quantitative prediction of the current through a molecule is still a great challenge 
despite the significant progress achieved in recent years [16-19]. For a more qualitative 
understanding of the mechanisms involved in molecular transport, it is thus advantageous 
to employ for the molecule a rather generic tight-binding model [8, 9, 20-25]. Then a 
flexible method for the computation of transport properties is provided by master equations 
of the Bloch-Redfield type which allow one to include electron-electron and electron-phonon 
interactions, as well as time-dependent fields [9,26]. Similar methods have also been used for 
describing incoherent transport [27,28]. 

Here, we explore the role of slow fluctuations or static disorder for molecular 
conductance. Thereby we will assume that the relevant environmental fluctuations are so slow 
that they can be described as static disorder which defines an ensemble of wire Hamiltonians. 
Then a natural quantity of interest is the corresponding distribution of stationary currents. 
A setup for which this current distribution is also directly relevant is an array of molecular 
junctions that conduct in parallel. We employ a tight-binding model for the molecule and 
treat it with the Floquet master equation formalism derived in Ref. [26] which we review 
briefly in Section 2. In Section 3, we present results for a static model with a large voltage 
bias, while in Section 4, we investigate pumping effects caused by an interplay of ac driving 
fields and disorder. The analytical derivation of the current distribution for a wire with two 
sites is deferred to the Appendix. 
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Figure 1. Bridged molecular wire model consisting of N = 5 sites with internal tunnelling 
matrix elements A and effective wire-lead coupling strengths r L /R. 



2. Wire-lead model and master equation 

The system of the driven molecular wire, the leads, and the coupling between the molecule 
and the leads, as sketched in Fig. 1, is described by the Hamiltonian 

•ffl (t) = J%wire(t ) ~\~ ^feads "I - ^wire— leads- (1) 

The wire is modelled by N tight-binding orbitals \n), n= 1, . . . ,N, such that 

N-l TJ 
^wrre = + ^n)c\c n - A £ [c\^C n + cjc+l ) + "1), (2) 

with the tunnel matrix element A and the capacitive energy U . Each onsite energy E n (t) + £, n 
contains a random contribution that subsumes the influence of environmental fluctuations. 
We assume that these fluctuations are Gaussian distributed and so slow that we can treat them 
as static disorder. Thus, the probability that the onsite energy of orbital n lies in an interval of 
size dt, around E n (t) + t, n reads 

where the variance a 2 is assumed to be position-independent. This implies that the energy 
fluctuations are spatially uncorrelated, such that (£ n £ n r) = <J 2 8 nn i. The onsite energies 
E n (t) = E n + Ax n cos(£lt) are modulated by a harmonically time-dependent dipole force, 
where A denotes the electrical field amplitude multiplied by the electron charge and the 
distance between neighbouring sites, with x n = j(N + 1 — 2n) the scaled position of site \n). 
Our goal will be to compute for many realizations of the wire Hamiltonian the resulting dc 
current which provides the current distribution P(I). The last term in Eq. (2) captures the 
electron-electron interaction within a capacitor model and the operator JV = c\c n describes 
the number of excess electrons residing on the molecule. Below we shall assume that U is so 
large that only states with zero or one excess electron play a role. The first and the last site of 
the molecule, 1 1) and \N), couple via the tunnelling Hamiltonian 

^ire-leads = £(^L 9 4 9 C 1 + V ^q C Rq C N) + H ' C ' ( 4 ) 
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to the respective lead. The operator c^ q (c Rq ) creates an electron in the left (right) lead in state 
\hq) which is orthogonal to all wire states. The influence of the tunnelling Hamiltonian is 
fully characterised by the spectral density T?{e) = IftY^q | V^^ | 2 5 (e — £ q ). If the lead states are 
dense and located at the centre of the conduction band, the spectral densities can be replaced 
by a constant, i.e. we assume r^(e) = T for both leads. 
The leads are modelled as ideal Fermi gases 

Pleads = £ (e ?L cl q C hq + £ qR cj^CjJ , (5) 

which are initially at thermal equilibrium with the chemical potential jU L / R and, thus, are 
described by the density operator pi ea ds,eq exp[— (J£f e ads _ Ml- 7 ^ — Mr-^R/'Ab^]) where 
,yp£ = Y, q c \fqi is the electron number operator for lead i = L, R. Since a typical metal 
screens all electric fields with a frequency below the plasma frequency, we assume that the 
bulk properties of the leads are not affected by the laser irradiation. 

2.1. Perturbation theory 

The derivation of a master equation starts from the Liouville-von Neumann equation ihp (t) = 
[Jf?(t),p(t)] for the total density operator p(t) for which one obtains by standard techniques 
the approximate equation of motion [9, 29-32] 

i 1 f°° ~ 

PM = -r[#wire(0 + #leads,p(0] ~ TT / dT [^wire-leads, [^wire-leads (* - T, t), p{t)}} . (6) 

n h Jo 

Here the first term corresponds to the coherent dynamics of both the wire electrons and the 
lead electrons, while the second term describes resonant electron tunnelling between the leads 
and the wire. The tilde denotes operators in the interaction picture with respect to the molecule 
and the lead Hamiltonian without the molecule-lead coupling, X(t,t') = ul(t,t')XUQ(t,t'), 
where Uq is the propagator without the coupling. The net (incoming minus outgoing) 
electrical current through the left contact is given by minus the time-derivative of the electron 
number in the left lead multiplied by the electron charge —e. From Eq. (6) follows for the 
current in the wide-band limit the expression 

I L (t)=etr[p(t),Ji] = -e^Re fdr / dee 1 ^ 

nh Jo J 

x { (cjcj (t, t - T))/ L (e) - (c i (t, t - T)cl)/ L (e) } , (7) 

where fy is the Fermi function of the respective lead and fg = 1 — fy. Furthermore, (• ■ ■) = 
tXwire p w j re • • • denotes the expectation value with respect to the wire density operator. We 
emphasise that the expectation values in Eq. (7) depend directly on the dynamics of the 
isolated wire and are thus influenced by the driving. 

2.2. Floquet theory 

An important feature of the current formula (7) is its dependence on solely the wire operators 
while all lead operators have been eliminated. Therefore it is sufficient to consider the reduced 
density operator p w i re = tri ea( js p of the wire electrons. The effort necessary to calculate p w i re 
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can be reduced significantly by exploiting the fact that the master equation (6) inherited from 
the total Hamiltonian Jj?(t) a periodic time-dependence, which opens the way for a Floquet 
treatment. 

One possibility for such a treatment is to use the Floquet states of the central system, i.e. 
the driven wire, as a basis. Thereby we also use the fact that in the wire Hamiltonian (2), 
the single-particle contribution commutes with the interaction term and, thus, these two 
Hamiltonians possess a complete set of common eigenstates. In analogy to the quasimomenta 
in Bloch theory for spatially periodic potentials, the quasienergies e a come in classes e a _k = 
e a +khQ., k E Z, of which all members represent the same physical solution of the Schrodinger 
equation. Thus we can restrict ourselves to states within one Brillouin zone like for example 
< e a < hCl. 

For the computation of the current it is convenient to have an explicit expression for 
the interaction picture representation of the wire operators. It can be obtained from the 
(fermionic) Floquet creation and annihilation operators defined via the transformation c a (t) = 
'L n (<Pa(t)\n)c n , which reads in the interaction picture c a (t,t') = e -i( £ a+u.yK, iK )(t-t')/h Ca ^f^ 
with the important feature that the time difference t — t' enters only via the exponential 
prefactor [9]. 



2.3. Master equation and current formula 

In the following, we assume the interaction U to be the dominant energy scale in the system, 
therefore we can restrict the wire Hilbert space to the N + 1 dimensional subspace of states 
with zero and one electron, such that a basis for the decomposition of the reduced operator is 
{|0),Ccj(/) |0)}, where |0) denotes the wire state in the absence of excess electrons. Moreover, 
it can be shown [26] that at large times, the density operator of the wire becomes diagonal in 
the electron number ,jV . Therefore a proper ansatz reads 

Pwire(0 = |0)P00(0(0| + £4|0)p^(0(0|^. (8) 

Note that we keep terms with a ^ j8, which means that we work beyond a rotating-wave 
approximation. 

Following our evaluation of the master equation [26], we arrive at a set of N 2 coupled 
equations of motion for p a p (t) which in Fourier representation read 

K e a- e p- kn &)Pap,k= -b £ (^a,F+F|l)(l|<P J 8l+F)Poo,F(/L(eal'+F)+/L(e j 8i + F)) 

z k',k" 

H (9a,k'+k"\^){M9a',k+k")Pa'p,k'f^( £ a',k+k") 
Z a',k',k" 

£ (Vp>,k>+k>>\l)(MVp,k+k'')Pal3>,k>fU £ P>,k>+k>>) 
P',k',k" 

+ same terms with the replacement 1 , L — > iV, R. (9) 
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In an analogous manner we obtain for the dc current the expression 
h = ^ Re £ ( £ ( (pp # +k \l)(l\(pa,k)Pap ,k'f^a,k) 

~ £(<Pa,*'+fcl 1) (1 1 (pa.k)POO.k'fU£a,k)) ■ (10) 
k' ' 

The results of this section allow us to numerically compute the dc current through a 
driven conductor as well as studying the undriven limit. The current distributions discussed 
below are obtained by computing the dc current for typically 10 4 realizations of the wire 
Hamiltonian (2). Then these values are taken for a histogram with 150 bins which finally will 
be scaled such that we obtain a normalised probability density. 

3. Electron transport with slowly fluctuating energies 

We first address an undriven wire in the configuration sketched in Fig. 1 where the distribution 
of all wire levels is centred at energy E„ = 0. The transport voltage is so large that 
all eigenenergies lie well within the voltage window and, consequently, the transport is 
unidirectional. Then in the absence of onsite energy fluctuations (o = 0), the current can be 
evaluated analytically within a rotating-wave approximation and reads 7 max = eT/h(N+ 1), 
i.e. it decays with increasing wire length [33]. The index "max" refers to the fact that any 
modification of the onsite energies can only reduce the current — which is confirmed by our 
simulations. The physical reason for this is that for equal onsite energies, solely the kinetic 
energy determines the eigenstates which, consequently, are delocalised. Different onsite 
energies, by contrast, tend to "localise" the eigenstates. Thus in the limit of small disorder, 
the current distribution P(I) is expect to possess a clear peak at / = 7 max and some minor 
contribution for lower values of /. 

Figure 2 shows the simulated current distributions for two different variances. For a 
small variance (panel a), the distributions for short wires show the expected behaviour. With 
an increasing wire length, the peak at / = 7 max disappears and is eventually replaced by an 
apparently parabolic distribution. This length dependence can be understood in the following 
way: For a short wire, the probability that a level is out of resonance is rather low and, thus, 
most realizations of the wire Hamiltonian will allow resonant inter-site tunnelling. With an 
increasing number of levels, however, the probability for having at least one misaligned level 
increases and a current significantly smaller than 7 max becomes more likely. The precise values 
will depend on the details and, consequently, we expect a broad distribution. This means that 
whenever a large number of levels plays a role, the transport through a molecule is extremely 
sensitive to even small disorder induced by environmental fluctuations. 

With a larger variance, this scenario becomes even more pronounced as can be seen in 
Fig. 2b: Then the peak at 7 max is rather small and noticeable only for 2 and 3 sites. The 
most likely realization is a completely disordered wire with an accordingly low conductance. 
For N > 3, the distributions even possesse a significant peak at / = which corresponds 
to isolating behaviour. A closer inspection of the numerical data reveals that the crossover 
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Figure 2. Current distribution for a channel with N sites in the limit of a large bias voltage. 
The standard deviation of the onsite energies is a = 0.5A (a) and a — 2A (b), while the wire- 
lead coupling is r = 0.1 A. The distributions have been obtained by computing the current 
for 1.5 x 10 4 realizations of the wire Hamiltonian. The black dotted lines mark the analytical 
results for N = 2 sites. 

between conducting and isolating behaviour occurs when the effective disorder y/No exceeds 
the tunnel matrix element A. 

Interestingly enough, for N = 2, 3 the distribution turns out to be even non-monotonic, 
which means that most one most likely finds either a current close to the theoretical maximum 
or a significantly smaller lower value. The non-monotonic behaviour for N = 2 can also 
be found analytically. The derivation of the corresponding current distribution (A. 3) can be 
found in the Appendix. The excellent agreement of this analytical solution and the simulated 
distributions emphasis that the simulation with approximately 10 4 realisations ensures good 
convergence. 
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Figure 3. Current distribution for an AC driven wire with N = 2 sites for various bias voltages. 
The fluctuations of the onsite energies are characterised by the standard deviation a = 0.5A, 
the driving frequency and amplitude are A = A and £2 = 2A/h, respectively. All the other 
parameters are as in Fig. 2. 

4. AC-driven disordered junctions 

In order to investigate the influence of an AC driving, we employ the same model as 
above, but now with an additional dipole driving modelled by time-dependent onsite energies 
E n (t) = Ax n cos(Q.t) as discussed in Sect. 2. The driving frequency Q. = 2A/h is chosen 
such that it matches the average splitting of the wire energies, while the amplitude A = A 
corresponds to intermediately strong driving. The solid line in Fig. 3 shows the current 
distribution in the absence of a voltage bias, V = 0. The reflection symmetry of the ensemble 
relates to the symmetric shape of the distribution, which implies that the current vanishes 
in the ensemble average. An individual realization of the wire, however, generally does not 
possess reflection symmetry because the random energy shifts are spatially uncorrected. This 
asymmetry in combination with the non-adiabatic driving induces a coherent ratchet current, 
i.e. a dc current even in the absence of any net voltage bias. In the present case, the ratchet 
current is of the order of 10-20 percent of the current observed above in the large bias limit. 
This order of magnitude is typical when the driving frequency or a multiple of the driving 
frequency lies close to an internal resonance, while the intensity is moderate [26]. In addition 
to the broad distribution of ratchet currents, P(I) exhibits a peak at / = 0. This stems from 
realizations for which the driving is well out of resonance. 

For a bias voltage V > 0, the ensemble no longer possesses reflection symmetry and 
the current distribution is shifted towards positive values (see Fig. 3). For sufficiently 
small voltages V < A/e, non-adiabatic pumping against the voltage bias is still possible. 
Rather surprisingly, the peak at zero current remains. It now corresponds to realizations for 
which on the one hand, the driving is off-resonant while on the other hand, both levels lie 
outside the voltage window. With an increasing bias voltage, the second condition is less 
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frequently fulfilled, and eventually the distribution assumes a form similar to that found for 
a large voltage in the absence of driving. Already for V ~ 4-A/e, the distribution is hardly 
distinguishable from the one shown in Fig. 2a for a wire with N = 2 sites in the absence of 
driving. 

5. Conclusions 

We have investigated the current through a molecular wire with disordered onsite energies. 
Such a disorder can stem from the interaction with slow fluctuations of background charges 
in the substrate. In particular, we computed the resulting current distribution for two typical 
cases, namely an "open transport channel" and a driven molecular wire for which random 
energy shifts break reflection symmetry and, thus, the driving can induce a ratchet current. 

The open transport channel is characterised by tight-binding levels with equal onsite 
energies, such that any misalignment stems from the disorder. Its main consequence is that 
as soon as the standard deviation of the onsite energies exceeds the tunnel matrix elements, 
the current distribution no longer peaks only at a finite value, but also at zero. For longer 
wires, only the peak at zero current remains. This isolating behaviour resembles Anderson 
localisation which is found for electrons in a one-dimensional disordered lattice [34]. Note 
however, that we here considered short wires far from the scaling limit in which Anderson 
localisation is usually studied. 

Since the random energy shifts break reflection symmetry, driving the molecular wire 
with a laser field induces ratchet currents for which we found a relatively broad distribution. 
If the driving frequency is far from any molecular excitation energy, the ratchet current will 
be rather small, and we indeed found that this is the case for very many wire realizations. 
It has the consequence that the corresponding distribution possesses a spike at zero current. 
This means that non-adiabatic pumping of electrons against a voltage bias becomes generally 
impossible whenever the relevant wire energy levels lie well within the voltage window. 

In conclusion, our results reveal that slow fluctuations or a static disorder can have a 
significant effect on molecular conduction. In various cases, the current distribution emerges 
to be rather flat, which means that even the magnitude of the current depends sensitively on 
environmental influences. 
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Appendix A. Analytical solution for two levels 

A wire model with N = 2 sites represents an analytically solvable case for which one 
observes a non-monotonic current distribution and which can serve as test case for numerical 
implementations. Here we consider a two-level system with on-site energies £1,2, i.e. with a 
bias 2r\ = t,\ — £2. Since the random energy shifts t, n are Gaussian distributed with variance 
a 2 , the bias 2r\ is also Gaussian distributed but with variance 2c 2 , i.e. its distribution function 
reads w(r\) = exp(-7] 2 /(7 2 ) /V no 2 . 

For the computation of the current, we restrict ourselves to the limit of a large transport 
voltage such that both eigenenergies of the two-level system lie within the voltage window. 
Then, the Fermi functions of the left and the right lead effectively become /l = 1 and 
/r = 0. In this case, transport can be described within rotating-wave approximation (RWA) 
which practically means that the reduced density operator of the wire is diagonal in energy 
representation [33]. Within RWA thus follows from the master equation (9) the occupation 
probability p aa = w l a /w 2 a and, thus, poo = 1 — La^/^cc The coefficients w n a = |(0 a |n)| 2 
denote the overlap between the eigenstate \<j) a ) and the localised state \n) (Note that in the 
undriven case, all non- vanishing contributions have sideband index k = 0, such that here the 
sideband index k can be omitted). Inserting this solution into the current formula (10), we 
obtain I = eF/h(l + E a w l a /w 2 a ) . 

The remaining task is now to diagonalise the single-particle Hamiltonian which provides 
the coefficients w n a . For bias 2r\ and tunnelling matrix element A, the Hamiltonian in pseudo- 
spin notation reads H — rio z + Ao x and possesses the eigenenergies ±8 = ±(r] 2 + A 2 ) 1 / 2 . 
The corresponding eigenvectors § a are proportional to ( 8 + r\ , A) and (5 — 17, A), respectively, 
such that w l a /w 2 a = (8 ± r/) 2 /A 2 . Then we obtain for the current the expression 

w \ ^ 1 ^max , . -. , 

W "l3+4T] 2 /A 2 ~ l+47] 2 /3A 2 ' { } 

which assumes its maximum 7 max = eT/3h in the unbiased limit r\ = 0. 
The probability distribution for the current relates to w(r\) via 

drii 



P(I) = Y,Hli) 



dl 



(A.2) 



where the summation considers all values of r\ that fulfil the condition / = I(f]). After some 
straightforward algebra, we obtain by evaluating expression (A.2) the current distribution 



/ 3A 2 / max // 2 / 3A 2 \ 

p{1) = Vw ^r=773T exp (-4^ (w// - !) )' (A - 3) 



4 ^Vw/-i 

which is defined and normalised on the interval [0, 7) 
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